Gene_clustering

Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls
3: In readChar(file, size, TRUE) : truncating string with embedded nuls

library(factoextra)
library(COTAN)
library(ggrepel)
library(Rtsne)
library(plotly)
library(tidyverse)
library(htmlwidgets)
library(MASS)
library(dendextend)
library(grid)
library(ggpubr)

To demostrate how to do the gene clustering usign COTAN we begin importing the COTAN object that stores all elaborated data and, in this case, regarding a mouse embrionic cortex dataset (developmental stage E17.5).

input_dir = "Data/"
layers = list("L1"=c("Reln","Lhx5"), "L2/3"=c("Satb2","Cux1"), "L4"=c("Rorb","Sox5") , "L5/6"=c("Bcl11b","Fezf2") , "Prog"=c("Vim","Hes1"))
#objE17 = readRDS(file = paste(input_dir,"E17.5_cortex.cotan.RDS", sep = ""))
objE17 = readRDS(file = paste(input_dir,"E17_cortex_cl2.cotan.RDS", sep = ""))

g.space = get.gene.coexpression.space(objE17, 
                                      n.genes.for.marker = 25,
                                      primary.markers = unlist(layers))

[1] "calculating gene coexpression space: output tanh of reduced coex matrix"
     L11      L12    L2/31    L2/32      L41      L42    L5/61    L5/62    Prog1    Prog2 
  "Reln"   "Lhx5"  "Satb2"   "Cux1"   "Rorb"   "Sox5" "Bcl11b"  "Fezf2"    "Vim"   "Hes1" 
[1] "Get p-values on a set of genes on columns genome wide on rows"
[1] "Using function S"
[1] "function to generate S "
[1] "Secondary markers:181"
[1] "function to generate S "
[1] "Columns (V set) number: 181 Rows (U set) number: 1236"

g.space = as.data.frame(as.matrix(g.space))

coex.pca.genes <- prcomp(t(g.space),
                 center = TRUE,
                 scale. = F) 

fviz_eig(coex.pca.genes, addlabels=TRUE,ncp = 10)

#fviz_eig(coex.pca.genes, choice = "eigenvalue", addlabels=TRUE)

Hierarchical clustering

par_pca = pca_1[colnames(g.space)[colnames(g.space) %in% rownames(pca_1)],]
#plot N 1
p1 <- ggparagraph(text = paste0(rownames(par_pca[par_pca$hclust == par_pca["Fezf2","hclust"],]), collapse = ", "), 
                  face = "italic", 
                  size =10, 
                  color = "#3C5488FF")

#plot N 2
p2 = ggparagraph(text = paste0(rownames(par_pca[par_pca$hclust == unique(par_pca$hclust)[!unique(par_pca$hclust) %in% unique(par_pca[unlist(layers),"hclust"])][1],]), collapse = ", "), 
                 face = "italic", 
                 size =10, 
                 color = "gray")

#plot N 3
p3 = ggparagraph(text = paste0(rownames(par_pca[par_pca$hclust == unique(par_pca$hclust)[!unique(par_pca$hclust) %in% unique(par_pca[unlist(layers),"hclust"])][2],]), collapse = ", "), 
                 face = "italic", 
                 size =10, 
                 color = "gray")
#plot N 4
p4 = ggparagraph(text = paste0(rownames(par_pca[par_pca$hclust == par_pca["Reln","hclust"],]), collapse = ", "), 
                 face = "italic", 
                 size =10, 
                 color = "#E64B35FF")

#plot N 5
p5 = ggparagraph(text = paste0(rownames(par_pca[par_pca$hclust == par_pca["Cux1","hclust"],]), collapse = ", "), 
                 face = "italic", 
                 size =10, 
                 color = par_pca["Cux1","colors"])
#plot N 6
p6 = ggparagraph(text = paste0(rownames(par_pca[par_pca$hclust == par_pca["Rorb","hclust"],]), collapse = ", "), 
                 face = "italic", 
                 size =10, 
                 color = par_pca["Rorb","colors"])
#plot N 7
p7 = ggparagraph(text = paste0(rownames(par_pca[par_pca$hclust == par_pca["Vim","hclust"],]), collapse = ", "), 
                 face = "italic", 
                 size =10, 
                 color = par_pca["Vim","colors"])

w = ggparagraph(text = " ", 
                 face = "italic", 
                 size =10, 
                 color = "white")

pp =ggarrange(p3,p5,p1,p2,p4,p6,p7,w,
          ncol = 1, nrow = 8,
          heights = c(0.1,0.15,0.23, 0.1, 0.2, 0.2, 0.1, 0.35))

    
lay <- rbind(c(1,NA),
             c(1,2.5),
             c(1,2.5),
             c(1,2.5),
             c(1,2.5),
             c(1,2.5),
             c(1,NA))

gridExtra::grid.arrange(tree, pp, layout_matrix = lay)

or just with primary markers

# some more genes as landmarks
controls =list("genes related to L5/6"=c("Foxp2","Tbr1"), "genes related to L2/3"=c("Mef2c"), "genes related to Prog"=c("Nes","Sox2") , "genes related to L1"=c() , "genes related to L4"=c()) 
dend %>%
dendextend::set("labels", ifelse(labels(dend) %in% rownames(pca_1)[rownames(pca_1) %in% c(unlist(layers),unlist(controls))], labels(dend), "")) %>%
  dendextend::set("branches_k_color", value = c("gray80","#4DBBD5FF","#91D1C2FF" ,"gray80","#F39B7FFF","#E64B35FF","#3C5488FF"), k = 7) %>%
plot(horiz=F, axes=T,ylim = c(0,100))

Now we can plot the PCA

# dataframe to be able to label only primary markers and control genes
textdf <- pca_1[rownames(pca_1) %in% c(unlist(layers),unlist(controls)) , ]

 for (m in c(1:length(controls))) {
  for (g in controls[[m]]) {
    if(g %in% rownames(textdf)){
      textdf[g,"highlight"] = names(controls[m])
    } 
  }
}

# deciding the colors
mycolours <- c("genes related to L5/6" = "#3C5488FF","genes related to L2/3"="#F39B7FFF","genes related to Prog"="#4DBBD5FF","genes related to L1"="#E64B35FF","genes related to L4" = "#91D1C2FF", "not marked"="#B09C85FF")

# to assing correcly the cluster number and the color
mycolours2 = c("Reln","Satb2","Rorb","Bcl11b","Vim")
names(mycolours2) = unique(cut[unlist(layers)])

mycolours2[mycolours2 == "Reln"] = "#E64B35FF"
mycolours2[mycolours2 == "Satb2"] = "#F39B7FFF"
mycolours2[mycolours2 == "Rorb"] = "#91D1C2FF"
mycolours2[mycolours2 == "Bcl11b"] = "#3C5488FF"
mycolours2[mycolours2 == "Vim"] = "#4DBBD5FF"
color_to_add = unique(pca_1$hclust)[!unique(pca_1$hclust) %in% as.numeric(names(mycolours2))]
names(color_to_add) = unique(pca_1$hclust)[!unique(pca_1$hclust) %in% as.numeric(names(mycolours2))]
color_to_add[color_to_add %in% 
                 unique(pca_1$hclust)[!unique(pca_1$hclust) %in% as.numeric(names(mycolours2))]] = "#B09C85FF"
mycolours2 = c(mycolours2,color_to_add)

pca1 = ggplot(subset(pca_1,!hclust %in% unique(cut[unlist(layers)])  ), aes(x=PC1, y=PC2)) +  geom_point(alpha = 0.3,color = "#B09C85FF",size=1)

#pca2 = pca1 + geom_point(data = subset(pca_1, highlight != "not marked" ), aes(x=PC1, y=PC2, colour=hclust),size=2.5,alpha = 0.9) 
pca2 = pca1 + geom_point(data = subset(pca_1, hclust %in% unique(cut[unlist(layers)]) ), aes(x=PC1, y=PC2, colour=as.character(hclust)),size=1,alpha = 0.5) + 
 scale_color_manual( "Status", values = mycolours2)  +
  scale_fill_manual( "Status", values = mycolours2)  +
  xlab("") + ylab("") +
  geom_label_repel(data =textdf , aes(x = PC1, y = PC2, label = rownames(textdf),fill=as.character(hclust)),
                   label.size = NA, 
                   alpha = 0.5, 
                   direction = "both",
                   na.rm=TRUE,
                   seed = 1234) +
  geom_label_repel(data =textdf , aes(x = PC1, y = PC2, label = rownames(textdf)),
                   label.size = NA, 
                   segment.color = 'black',
                   segment.size = 0.5,
                   direction = "both",
                   alpha = 1, 
                   na.rm=TRUE,
                   fill = NA,
                   seed = 1234) +
  ggtitle("PCA") +
  theme_light(base_size=10) +
    theme(axis.text.x=element_blank(),plot.title = element_text(size=14, 
                                    face="italic", 
                                    color="#3C5488FF",
                                    hjust=0.01,
                                    lineheight=1.2,margin = margin(t = 5, b = -15)),
        axis.text.y=element_blank(),
        legend.position = "none")  # titl)

pca2 #+ geom_encircle(data = pca_1, aes(group=hclust)) 

t-SNE code and plot

# run the t-SNE
cl.genes.tsne = Rtsne(g.space ,initial_dims = 100, dims = 2, perplexity=30,eta = 200, verbose=F, max_iter = 3000,theta=0.4,num_threads = 10,pca_center = T, pca_scale = FALSE, normalize = T )

d_tsne_1 = as.data.frame(cl.genes.tsne$Y)
rownames(d_tsne_1) = rownames(g.space)

d_tsne_1 = d_tsne_1[order.dendrogram(dend),]

# save the cluster numebr inside a dataframe with the t-SNE information
d_tsne_1$hclust = cut

d_tsne_1$names = rownames(d_tsne_1)

# as before to label only some genes
textdf <- d_tsne_1[rownames(d_tsne_1) %in% c(unlist(layers),unlist(controls)),]

for (m in c(1:length(controls))) {
  for (g in controls[[m]]) {
    if(g %in% rownames(textdf)){
      textdf[g,"highlight"] = names(controls[m])
    } 
  }
}


 p1 = ggplot(subset(d_tsne_1,!hclust %in% unique(cut[unlist(layers)])), aes(x=V1, y=V2)) +  geom_point(alpha = 0.3, color = "#B09C85FF", size=1)

p2 = p1 + geom_point(data = subset(d_tsne_1, hclust %in% unique(cut[unlist(layers)]) ), aes(x=V1, y=V2, colour=as.character(hclust)),size=1,alpha = 0.5) +
    scale_color_manual("Status", values = mycolours2)  +
  scale_fill_manual("Status", values = mycolours2)  +
  xlab("") + ylab("")+
  geom_label_repel(data =textdf , aes(x = V1, y = V2, label = names,fill=as.character(hclust)),
                   label.size = NA, 
                   alpha = 0.5, 
                   direction = "both",
                   na.rm=TRUE,
                   seed = 1234) +
  geom_label_repel(data =textdf , aes(x = V1, y = V2, label = names),
                   label.size = NA, 
                   segment.color = 'black',
                   segment.size = 0.5,
                   direction = "both",
                   alpha = 1, 
                   na.rm=TRUE,
                   fill = NA,
                   seed = 1234) +
  ggtitle("t-SNE") +
  theme_light(base_size=10) +
  theme(axis.text.x=element_blank(),plot.title = element_text(size=14, 
                                    face="italic", 
                                    color="#3C5488FF",
                                    hjust=0.01,
                                    lineheight=1.2,margin = margin(t = 5, b = -15)),
        axis.text.y=element_blank(),
        legend.position = "none")  # titl)
p2

Code to create an iteractive plot. This can be modified to be used with all the plots.

p = ggplot(d_tsne_1, aes(x=V1, y=V2, text= paste("gene: ",names))) +  
  geom_point(size=2, aes(colour=as.character(hclust)), alpha=0.8) +
  scale_color_manual("Status", values = mycolours2)  +
  xlab("") + ylab("") +
  ggtitle("t-SNE") +
  theme_light(base_size=10) +
  theme(axis.text.x=element_blank(),
        axis.text.y=element_blank())

ggplotly(p)

Multidimensional scaling (MDS) and plot

sessionInfo()

R version 4.0.4 (2021-02-15)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.5 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] ggpubr_0.4.0      dendextend_1.14.0 MASS_7.3-53.1     htmlwidgets_1.5.3 forcats_0.5.1     stringr_1.4.0    
 [7] dplyr_1.0.4       purrr_0.3.4       readr_1.4.0       tidyr_1.1.2       tibble_3.0.6      tidyverse_1.3.0  
[13] plotly_4.9.3      Rtsne_0.15        ggrepel_0.9.1     COTAN_0.1.0       factoextra_1.0.7  ggplot2_3.3.3    

loaded via a namespace (and not attached):
 [1] matrixStats_0.58.0   fs_1.5.0             lubridate_1.7.9.2    filelock_1.0.2       RColorBrewer_1.1-2  
 [6] httr_1.4.2           tools_4.0.4          backports_1.2.1      R6_2.5.0             DBI_1.1.1           
[11] lazyeval_0.2.2       BiocGenerics_0.36.0  colorspace_2.0-0     GetoptLong_1.0.5     withr_2.4.1         
[16] tidyselect_1.1.0     gridExtra_2.3        curl_4.3             compiler_4.0.4       cli_2.3.0           
[21] rvest_0.3.6          Cairo_1.5-12.2       basilisk.utils_1.2.2 xml2_1.3.2           labeling_0.4.2      
[26] scales_1.1.1         rappdirs_0.3.3       digest_0.6.27        foreign_0.8-81       rmarkdown_2.7       
[31] rio_0.5.16           basilisk_1.2.1       pkgconfig_2.0.3      htmltools_0.5.1.1    dbplyr_2.1.0        
[36] rlang_0.4.10         GlobalOptions_0.1.2  readxl_1.3.1         rstudioapi_0.13      gridGraphics_0.5-1  
[41] farver_2.0.3         shape_1.4.5          generics_0.1.0       jsonlite_1.7.2       crosstalk_1.1.1     
[46] zip_2.1.1            car_3.0-10           magrittr_2.0.1       Matrix_1.3-2         Rcpp_1.0.6          
[51] munsell_0.5.0        S4Vectors_0.28.1     abind_1.4-5          reticulate_1.18      viridis_0.5.1       
[56] lifecycle_0.2.0      stringi_1.5.3        yaml_2.2.1           carData_3.0-4        parallel_4.0.4      
[61] crayon_1.4.0         lattice_0.20-41      haven_2.3.1          cowplot_1.1.1        circlize_0.4.12     
[66] hms_1.0.0            knitr_1.31           ComplexHeatmap_2.6.2 pillar_1.4.7         rjson_0.2.20        
[71] ggsignif_0.6.1       stats4_4.0.4         reprex_1.0.0         glue_1.4.2           evaluate_0.14       
[76] data.table_1.13.6    modelr_0.1.8         png_0.1-7            vctrs_0.3.6          cellranger_1.1.0    
[81] gtable_0.3.0         clue_0.3-58          assertthat_0.2.1     openxlsx_4.2.3       xfun_0.20           
[86] broom_0.7.5          rstatix_0.7.0        viridisLite_0.3.0    IRanges_2.24.1       cluster_2.1.1       
[91] ellipsis_0.3.1